function vmd_Espectra(fbase, fnum)
%
%   vmd_Espectra(fbase, fnum)
%
%   Reads file with density distribution and file with x-velocity.
%   Creates bunch of images.
%
%   Input:
%        fbase   base of file number (runname)
%        fnum    file number
%
%   Domain information and colormap limits are edited in the file.



%---------------------------------------------------------

% system size

jx=48; jy=1; jz=16;

% physics paramters

Lx=450; Ly=Lx; kz=0.34; nz=1; 


%---------------------------------------------------------


dx = Lx/jx;
dy = Ly/jy;
dz = 2*pi*nz/(jz*kz);

dkz = kz/nz;
dkx = 2*pi / Lx;
dky = 2*pi / Ly;

ixc = jx/2+1;
iyc = min(jy/2+1, jy);

%---------------------------------------------------------

%-- reading data --

  fn = num2str(fnum,'%04d');

  fname = [fbase, '.Ex.', fn];
  fid = fopen(fname, 'rb');
  f = fread(fid, jx*jy*jz, 'double');
  fclose(fid);
  Ex = reshape(f, jz,jy,jx);

  fname = [fbase, '.Ey.', fn];
  fid = fopen(fname, 'rb');
  f = fread(fid, jx*jy*jz, 'double');
  fclose(fid);
  Ey = reshape(f, jz,jy,jx);

  fname = [fbase, '.Ez.', fn];
  fid = fopen(fname, 'rb');
  f = fread(fid, jx*jy*jz, 'double');
  fclose(fid);
  Ez = reshape(f, jz,jy,jx);

  clear f

%-- find maximum of |E|^2 ---

  Esq = Ex.*Ex + Ey.*Ey + Ez.*Ez;

  Esq_avg = sum(sum(sum(Esq))) / (jx*jy*jz);

  [Esq_max, zmax] = max(Esq); 
  [Esq_max, ymax] = max(Esq_max); 
  [Esq_max, xmax] = max(Esq_max); 

  ix=xmax(1,1);
  if (jy == 1)
    iy = 1;
  else
    iy=ymax(1,1,ix);
  end
  iz=zmax(1,iy,ix);

  fprintf('Esq_avg = %12.6e\n',  Esq_avg);
  fprintf('Esq_max = %12.6e\n',  Esq_max);

  fprintf('Emax    = %12.6e,  at ix = %i, iy = %i, iz = %i\n',  ...
	  sqrt(Esq_max), ix, iy, iz);


%-- screen: spectrum of Ex(x, y=0, z=zmax) --

  E = Ex(iz,iy,:);
  E = reshape(E,1,jx);

  Ek = fft(E)/(jx);
  Eksq = Ek.*conj(Ek);

  figure(1, "Name", "Ex spectrum")
    plot ([0:jx/2-1]*dkx, Eksq(1:jx/2), '-o' )
  xlabel("k_x"); ylabel("|E_x(k_x)|^2")


%-- file: z-averaged spectra of Ex, Ey, and Ez along y=0 --

  Ek   = fftn( Ex(:,iyc,:) )/(jx*jz);
  Eksq = Ek.*conj(Ek);
  Exsq = sum(Eksq,1)/jz;

  Ek   = fftn( Ey(:,iyc,:) )/(jx*jz);
  Eksq = Ek.*conj(Ek);
  Eysq = sum(Eksq,1)/jz;

  Ek   = fftn( Ez(:,iyc,:) )/(jx*jz);
  Eksq = Ek.*conj(Ek);
  Ezsq = sum(Eksq,1)/jz;


  fname = [fbase, '.EfieldFTx.', fn, '.txt'];
  fid = fopen(fname, 'wt');

  fprintf(fid,'%% z-averaged x-spectra of E-field\n');
  fprintf(fid,'%% 1.kx   2.<|Ex(kx)|^2>   3.<|Ey(kx)|^2>   4.<|Ez(kx)|^2>\n')

  for i=1:jx/2
    k=(i-1)*dkx;
    fprintf(fid, "%14.6e %14.6e %14.6e %14.6e\n", k, Exsq(i), Eysq(i), Ezsq(i));
  end
  fclose(fid);



%-- screen: spectrum of Ez(x=0, y=0, z) --

  E = Ez(:,iyc,ixc);
  E = reshape(E,jz,1);

  Ek = fft(E)/(jz);
  Eksq = Ek.*conj(Ek);

  figure(2, "Name", "Ez spectrum")
  plot ([0:jz/2-1]*dkz, Eksq(1:jz/2), '-o' )
  xlabel("k_z"); ylabel("|E_z(k_z)|^2")


%-- file: spectrum of Ez(x=0, y=0, z) --


  fname = [fbase, '.EfieldFTz.', fn, '.txt'];
  fid = fopen(fname, 'wt');

  fprintf(fid,'%% z-spectrum of Ez along center line\n');
  fprintf(fid,'%% 1.kz   2.<|Ez(kz)|^2>\n')

  for i=1:jz/2
    k=(i-1)*dkz;
    fprintf(fid, "%14.6e %14.6e\n", k, Eksq(i));
  end
  fclose(fid);



%---------------------------------------------------------

end
